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We present results for the matrix elements of the additional AS — 2 operators that appear in 
models of physics beyond the Standard Model (BSM), expressed in terms of four BSM B-parameters. 

Combined with experimental results for AMk and ck, these constrain the parameters of BSM mod¬ 
els. We use improved staggered fermions, with valence HYP-smeared quarks and Nf = 2-1-1 flavors 
of “asqtad” sea quarks. The configurations have been generated by the MILC collaboration. The 
matching between lattice and continuum four-fermion operators and bilinears is done perturbatively 
at one-loop order. We use three lattice spacings for the continuum extrapolation: a ~ 0.09, 0.06 
and 0.045 fm. Valence light-quark masses range down to « while the light sea-quark 

masses range down to ~ mg®^®/20. Compared to our previous published work, we have added four 
additional lattice ensembles, leading to better controlled extrapolations in the lattice spacing and 
sea-quark masses. We report final results for two renormalization scales, p = 2 GeV and 3 GeV, 
and compare them to those obtained by other collaborations. Agreement is found for two of the 
four BSM B-parameters (B 2 and Bf'^®®^). The other two (B 4 and B 5 ) differ significantly from those 
obtained using RI-MOM renormalization as an intermediate scheme, but are in agreement with re¬ 
cent preliminary results obtained by the RBC-UKQCD collaboration using RI-SMOM intermediate 
schemes. 
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I. INTRODUCTION 

Neutral kaon mixing and the associated indirect CP- 
violation have long provided an important window into 
physics at high energy scales. In the Standard Model 
(SM), for example, the measured CP-violating parame¬ 
ter ck is sensitive to scales up to the top-quark mass. To 
determine whether the measured value is consistent with 
the SM, however, requires knowledge of the hadronic ma¬ 
trix element parametrized by the kaon B-parameter, Bk- 
Recently, lattice QCD calculations have matured to the 
point that such matrix elements can be determined from 
first principles with percent-level accuracy}^ Specifically, 
results for Bk from Refs. [m are such that the average 
has an error of ^ 1.3% [T]. This is accurate enough to 
provide strong constraints on SM parameters (see, e.g.. 
Refs, [ini [TT]1. Ultimately, lattice calculations will also 
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^ For a recent review of such quantities and their associated errors, 
see Ref. [I]. 


be able to use the RTl — Ks mass difference, AMk, to 
test the SM [T^ . 

Physics beyond the SM (BSM) will, in general, con¬ 
tribute to flavor changing neutral processes such as kaon 
mixing. Indeed, unless there is some cancellation akin 
to the GIM mechanism, rough estimates show that the 
scale of new physics must be > 10® TeV in order to 
avoid overly large contributions to AMk and ek [13] ■ In 
fact, many BSM models have partial cancellations such 
that the scale of new physics is accessible at the Large 
Hadron Collider (LHC), but often such models are push¬ 
ing against the constraints from kaon mixing. If evidence 
for new physics is discovered at the LHC in the coming 
years, then, in order to sift through the available mod¬ 
els, it will be essential to turn the constraints from kaon 
mixing into precision tools. To do this it is necessary to 
calculate the hadronic matrix elements of the full basis 
of AS = 2 four-fermion operators that can appear. Il¬ 
lustrations of how these matrix elements constrain BSM 
models are given in Refs. [HHii]- 

In the SM, four-fermion operators in the effective 
AS = 2 Hamiltonian are composed of left-handed cur¬ 
rents. Generic BSM physics, by contrast, also includes 
heavy virtual particles coupling to right-handed quarks. 
Because of this, the single “left-left” AS = 2 four-fermion 
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operator is augmented by four additional operators. Our 
aim in the present work is to provide fully controlled 
results for the corresponding additional mixing matrix 
elements. 

Calculations of such matrix elements using lattice 
QCD have a fairly long history. Initial results were ob¬ 
tained starting in the late 1990s in the quenched approx¬ 
imation [TMSn] . Then, in 2012, first results with un¬ 
quenched light quarks were presented by the ETM [13] 
and RBC-UKQCD collaborations jUj. These calcula¬ 
tions used, respectively, twisted-mass and domain-wall 
lattice fermions. Both performed the matching of lattice 
and continuum operators using non-perturbative renor¬ 
malization (NPR) ^21 and the RI-MOM scheme. The 
results for all four BSM B-parameters were consistent 
between the two calculations. 

In 2013, we presented results from a first calcula¬ 
tion of the BSM B-parameters using improved staggered 
fermions and one-loop perturbative matching of lattice 
and continuum operators |53|. Our results disagreed sig¬ 
nificantly for two of the four B-parameters with those 
from Refs. Haul]. In 2014, we discovered a minor error 
in our analysis that changed our results by ^ 5%. We 
also extended the range of lattice ensembles studied, so 
that the continuum and chiral extrapolations were bet¬ 
ter controlled. Preliminary results correcting the analysis 
and incorporating the new ensembles were presented in 
Ref. [21|. The discrepancy with Refs. [nm] remained 
at about the 3cr level for two of the B-parameters. 

The purpose of the present paper is provide a detailed 
description of our calculation along with our final results. 
In fact, these results are very close to the preliminary 
numbers presented in Ref. |24] , but there are many details 
not provided in either Ref. [23] or [24] that we present 
here. A further motivation for this work is provided by 
the recent results from the ETM and RBC-UKQCD col¬ 
laborations, presented in Refs. |3] and at Lattice 2015 
[25l [26] , respectively. The former work (which extends 
the Nf = 2 simulations of Ref. m to N — 2 -f 1 -p 1) es¬ 
sentially confirms the earlier results of Ref. [T3] , and thus 
continues to disagree with our results. The latter calcu¬ 
lation, Ref. [IS] , presents an investigation of the origin of 
the discrepancies by repeating their computation with a 
second lattice spacing and performing the renormaliza¬ 
tion with various schemes, including RI-SMOM schemes 
with non-exceptional kinematics [57]. Although the dis¬ 
cretization artifacts are found to be larger than previ¬ 
ously anticipated, the most important effects come from 
the renormalization procedure. The preliminary results 
of RBC-UKQCD with the new SMOM schemes are in 
approximate agreement with those presented here |26j . 
Given this complicated and confusing situation, it is im¬ 
portant to have a clear description of the details of all 
the calculations. 

Our work relies on several auxiliary theoretical calcula¬ 
tions. For the chiral extrapolations we need results from 
SU(2) staggered chiral perturbation theory (SChPT), 
and these are provided in Ref. [55] . We also need to know 


how to set up the calculation using staggered fermions 
(i.e. dealing with the extra valence tastes) as well as 
one-loop matching factors. These results are provided in 
Refs. [551I5D]. Finally, we need to evolve matrix elements 
using the continuum renormalization group for AS” = 2 
operators. The required two-loop anomalous dimensions 
were calculated in Ref. m, and some additional techni¬ 
cal details are worked out in Ref. [50] . 

This paper is organized as follows. In Sec. [ll] we de¬ 
scribe the basis of AS = 2 four-quark operators that 
we use, and the corresponding B-parameters and gold- 
plated combinations. In Sec. m we describe the details 
of the lattice calculation. We next turn to the analy¬ 
sis. Sec. jlVj explains how we extrapolate valence quark 
masses to their physical values, while Sec. [V] describes 
the combined extrapolation to the continuum limit and 
to physical sea-quark masses. We present our final re¬ 
sults and error budget in Sec. [VTJ and compare these to 
the above-mentioned results that use other fermion dis¬ 
cretizations in Sec. m 

II. AS = 2 FOUR-QUARK OPERATORS AND 
BAG PARAMETERS 

We use the operator basis (Buras’s basis) of Ref. [5T] . 
in which the AS = 2 four-quark operators are 

Qi = [s“7m( 1 - 75)d“][s%M(l - 75)d''], 

Q2 = r(l-75)d1[s'’(l-75)d'’], 

Qs = - 75)rf“][s'’cr^!.(l - l5)d^] , (1) 

Q4 = r(l-75)d1[s''(l + 75)d''], 

Qd = [s“7/x(l - 75)d“][s%^.(l -k75)d^]. 

Here the operators have been written in Euclidean space, 
with a and b being color indices. Qi is the operator cor¬ 
responding to Bk, while <32,3,4,5 are the BSM operators^ 

The hadronic matrix elements of the AS = 2 four- 
quark operators can be parametrized by so-called kaon 
bag parameters (or B-parameters). These are conven¬ 
tionally defined by 

_ ( 2 ) 

|(Aro|s7^75d|0)(0|s7^75d|A:o) 

^ ^ _ (Ko|Q,-|Ko) 

' iV,(Ko|s75d|0)(0|s75d|Ko) ’ 

where j = 2 — 5, and 

(7V2, iVa, A4, As) = (5/3, 4, -2, 4/3) (4) 


^ This basis is complete in four dimensions aside from the need to 
add the parity conjugates of Qi — Q 3 . We do not consider these 
additional operators, however, since they have the same positive 
parity parts as Qi — Q 3 , and the matrix element we consider 
picks out the positive parity parts. 
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are factors arising in the vacuum saturation approxima¬ 
tion. In the following, we will often refer collectively 
to “the Bi”, and this will indicate all five of the B- 
parameters, i.e. the index i runs over i = K,2, 3,4, 5. 

In our lattice calculation, we find it more convenient 
to evaluate the B-parameters rather than the correspond¬ 
ing matrix elements, {Ko\Qi\Ko). This avoids the need 
to determine the overlap of our sources with the Kq and 
Kq states, reduces the dependence on the scale, a, since 
the B-parameters are dimensionless, cancels some of sta¬ 
tistical and systematic errors, and simplifies chiral ex¬ 
pansions, since the staggered chiral perturbation theory 
(SChPT) expressions are simpler pH] , 

We also make extensive use of “gold-plated” combi¬ 
nations of the B-parameters. These are combinations 
chosen to be free of chiral logarithms at next-to-leading 
order (NLO) in SU(2) chiral perturbation theory |28) : 


G 21 

G24 


B2 _ B2 

R ’ 23 D ’ 


B 2 • Bi, G 45 


Bb 


(5) 


In this paper, the subindex i of the Gi runs over i = 
21,23,24,45. 

As described below, it turns out that the combined ex¬ 
trapolation in a? and sea-quark masses is much better 
controlled for the Gi and Bk than for B 2 - 5 . Thus our fi¬ 
nal results for the BSM B-parameters are obtained using 
Bk and the Gi in the following way: 

B^ = Bk ■ G 21 , 
uG _ n G 21 

B3 - Bk - , 

^23 

„G _ G 24 ( 6 ) 

" Bk ■ G 21 ’ 

^ G 24 

^ Bk ■ G 21 • G 45 

The superscript G indicates that we use gold-plated com¬ 
binations to reconstruct the B-parameters. 


III. LATTICES AND MEASUREMENTS 

We use the MILC ensembles listed in TablelH These are 
generated with A/ = 2 - 1-1 flavors of staggered fermions 
using the “asqtad” fermion action. Details of the con¬ 
figuration generation are given in Ref. [32) . To convert 
our data to physical units, we use the values of ri/a 
obtained by the MILC collaboration. [32l |33], and set 
ri = 0.3117(22) fm, following Refs. [321l31]jf] We stress 


^ Some values of ri/a are updated compared to Ref. [32]; these are 
F2: 3.6987, F3: 3.7036, F4: 3.7086, F5:3.6993, F7: 3.7000, F9: 
3.6984, S4: 5.2825, S5: 5.2836 [33]. 


TABLE I. MILC ensembles used in our numerical study. Here 
“ens” represents the number of gauge configurations, “meas” 
is the number of measurements per configuration, and ID is 
a label, ami and arris are, respectively, the light and strange 
sea quark masses in lattice units. The values of a are nominal. 


a (fm) 

ami jams 

size 

ens X meas 

ID 

0.12 

0.03/0.05 

20^ X 64 

564 X 9 

Cl 

0.12 

0.02/0.05 

20^ X 64 

486 X 9 

C2 

0.12 

0.01/0.05 

20^ X 64 

671 X 9 

C3 

0.12 

0.01/0.05 

28^ X 64 

274 X 8 

C3-2 

0.12 

0.007/0.05 

20^ X 64 

651 X 10 

C4 

0.12 

0.005/0.05 

24^ X 64 

509 X 9 

C5 

0.09 

0.0062/0.031 

28^ X 96 

995 X 9 

FI 

0.09 

0.0031/0.031 

40^ X 96 

959 X 9 

F2 

0.09 

0.0093/0.031 

28^ X 96 

949 X 9 

F3 

0.09 

0.0124/0.031 

28^ X 96 

1995 X 9 

F4 

0.09 

0.00465/0.031 

32^ X 96 

651 X 9 

F5 

0.09 

0.0062/0.0186 

28^ X 96 

950 X 9 

F6 

0.09 

0.0031/0.0186 

40^ X 96 

701 X 9 

F7 

0.09 

0.00155/0.031 

64^ X 96 

790 X 9 

F9 

0.06 

0.0036/0.018 

48^ X 144 

749 X 9 

SI 

0.06 

0.0025/0.018 

56® X 144 

799 X 9 

S2 

0.06 

0.0072/0.018 

48® X 144 

593 X 9 

S3 

0.06 

0.0054/0.018 

48® X 144 

582 X 9 

S4 

0.06 

0.0018/0.018 

64® X 144 

572 X 9 

S5 

0.045 

0.0028/0.014 

64® X 192 

747 X 1 

U1 


that the values of a listed in the table are nominal. The 
actual values (determined from ri /a) differ slightly from 
the nominal values, and it is the former that we use in 
our analysis. In the following, we sometimes use MILC 
terminology and refer to the sets of ensembles with nom¬ 
inal lattice spacings of a = 0.12 fm, 0.09 fm, 0.06 fm and 
0.045 fm as coarse, fine, superfine and ultrafine lattices, 
respectively. 

Compared to the results presented in Ref. [21], the 
additional ensembles are F 6 , F7, F9 and S5. These ad¬ 
ditions significantly improve the reliability of the chiral 
extrapolations, as we now explain. On all ensembles ex¬ 
cept F 6 and F7, the strange sea quark masses lie close 
to, but not exactly at, the physical value. Adding in F 6 
and F7, which have lighter strange sea quarks, allows us 
to correct for the offset. Adding S5 ensures that on both 
fine and superfine lattices the average up/down sea quark 
mass, TO£, ranges down to « mP^^®/ 10 , so that the chiral 
extrapolation is relatively short. Finally, adding in F9 
provides us with a light sea quark mass, mi « m^^^^/20, 
that lies much closer to the physical value. 

We use HYP-smeared staggered fermions |3S] as va¬ 
lence quarks. Parameters for the HYP smearing are cho¬ 
sen to remove G(a^) taste-symmetry breaking at tree 
level [35]. We use 10 different valence quark masses on 
each lattice: 

Tl 

-rrix, my = to™™ x — with n = 1, 2, 3,..., 10 , 


( 7 ) 
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where m”°™ is the nominal strange quark mass given in 
Table im We have labeled the valence masses and Wy, 
the former corresponding to the valence d quark and the 
latter to the valence s quark. 

As explained in the next section, nix and ruy will be 
extrapolated to their physical values, and 

respectively. To determine these physical values on each 
ensemble use the same method as in Ref. [7]. First, the 
flavor non-singlet yy “pion” mass is extrapolated until is 
equals Mss,phys = 0.6858(40) GeV, which is the “phys¬ 
ical” value determined in Ref. m- This determines 
j.^phys_ is extrapolated (with my at its now- 

determined physical value) such that the xy “kaon” has 
a mass equal to that of the physical Kq. These extrapo¬ 
lations are done separately on each ensemble. For illus¬ 
tration, we show in Table |TT] the resulting physical values 
(as well as the valence masses we use in the simulations) 
for the ensembles having mi/m^ = 1/5. We see that our 
lightest valence quark masses are roughly twice 
while our heaviest lie somewhat below 


TABLE II. Physical values of valence quark masses on repre¬ 
sentative ensembles, in lattice units. For comparison, we also 
show the range of valence masses used in simulations. 


Ensemble 

phvs 

am)) ^ 


arrix and arriy 

C3 

0.00213(2) 

0.05204(5) 

0.005-0.05 

FI 

0.00146(2) 

0.03542(5) 

0.003-0.03 

SI 

0.00104(1) 

0.02372(3) 

0.0018-0.018 

U1 

0.00076(1) 

0.01693(3) 

0.0014-0.014 


We calculate the valence xx “pion” and xy “kaon” 
masses in standard fashion using the same wall sources 
as described below. The statistical errors on these results 
are very small. In Table [llT| we quote some representative 
values to indicate the range of pion masses in physical 
units. Note that (val, max) is the mass of the heav¬ 
iest pion that we use in our chiral extrapolation to the 
physical valence d quark. This extrapolation is discussed 
in the following section. We also include values for the 
lightest sea-quark pion for the fine, superfine and ultra- 
fine lattices, as well as for the coarse ensemble that we 
use to study finite-volume effects. 

We use essentially the same methodology for calculat¬ 
ing the BSM B-parameters as we employed in the cal¬ 
culation of Bk in Ref. [7]. Thus we give only a brief 
discussion here, while for Bx we refer to Ref. [7]. In 
terms of lattice operators, the BSM B-parameters are 

' iV,(Z°pi|zpOL-‘(t)|0)(0|^pOL-‘(t)|A:02) ’ 

where Q})**'* are lattice four-fermion operators and Op®'* 
is the taste ^5 pseudoscalar bilinear, zjk and zp are 
one-loop matching factors that convert lattice operators 
to their continuum counterparts, the latter defined in 
the MS scheme using naive dimensional regularization 


TABLE III. Valence and sea pion masses (in GeV) on repre¬ 
sentative ensembles, (val, min) and m^(val, max) are the 
minimum and maximum valence pion masses used in our va¬ 
lence chiral extrapolation. The values for these quantities on 
other coarse, fine and superfine ensembles are very similar to 
those on ensembles C3, F9 and S5, respectively. For the fine, 
superfine and ultrafine lattices, we show the sea quark pion 
mass on the ensemble with the smallest value of this quantity. 
For the coarse lattices, we pick the ensemble used to estimate 
finite-volume effects. 


Ensemble 

m,r(val, min) 

m.,r(val, max) 

rrin (sea) 

C3 

0.222 

0.430 

0.372 

F9 

0.206 

0.401 

0.174 

S5 

0.195 

0.379 

0.222 

U1 

0.206 

0.397 

0.316 


(NDR). We use the mean-field improved lattice operators 
defined in Refs. [7ni[3D]. The one-loop matching is quite 
involved as one must ensure that the continuum basis is 
extended to d = 4 — 2e dimensions using the same defini¬ 
tion of evanescent operators as in Ref. [ST]. The matching 
factors have been worked out and described in detail in 
Ref. ISO], building on the earlier work of Ref. |00j, and 
we do not repeat them here. They depend on the renor¬ 
malization scale y of the continuum operator and on Og. 
The latter is chosen to be in the MS scheme and is evalu¬ 
ated at the same scale /r. In our initial matching we take 
/i = I/a and then evolve the results in the continuum to a 
common renormalization scale. In the numerical evalua¬ 
tion of the matching coefficients we use four loop running 
to determine as{y), using as input as{Mz) = 0.118. 

To produce the kaons and antikaons, we place U(l)- 
noise wall sources on time slices ti and ^ 2 , with ^2 > ti- 
These produce taste ^5 kaons and anti-kaons having zero 
spatial momenta. The four-quark operators are placed 
between the sources at time t (i.e. ti < t < ^ 2 )- When 
t is far enough from the sources, so that excited state 
contamination is small, the three-point correlators should 
be independent of t, and can be fit to a constant. To 
determine the fit range, we use the two-point correlator 
from the wall-source to the taste ^5 axial current. From 
the effective mass plot for this correlator, we find the 
distance from the source, tp, for which the contamination 
from excited states becomes negligibly small. Then we 
fit from t = ti + tp to t = ti + tp = t 2 — tp — 1 (which is 
a symmetrical range since our operators extend over the 
two time slices t and t + 1). Our choices oi tp and tp 
are given in Table IV Note that we choose At = t 2 — ti 
to be less than half of the time extent of the lattice to 
avoid “around the world” contributions. Further details 
concerning sources and time ranges is given in Ref. [7]. 


The plateaus resulting from the above-described pro¬ 
cedure are reasonable. Examples are shown for the 
gold-plated combinations G23 and G45 in Figs. and 
Here we show cases with light valence quark masses 
{mxjmy = 1/10 with mi/ms = 1/5) for which the sta¬ 
tistical errors are larger. The fits to a constant are per- 
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TABLE IV. Choices for the wall source separation, At = t 2 — 
ti, and its ratio to the temporal length of the lattices, T, as 
well as the parameters determining the fitting range. 


lattice spacing 

At 

At/T 

th 

tR 

ti (fm) 

0.12 fm 

26 

0.41 

10 

15 

1.19 

0.09 fm 

40 

0.42 

14 

25 

1.18 

0.06 fm 

60 

0.42 

22 

37 

1.29 

0.045 fm 

80 

0.42 

26 

53 

1.14 


formed ignoring correlations between time slices (diago¬ 
nal approximation for the covariance matrix) in order to 
avoid instabilities due to the small eigenvalues of the co- 
variance matrix [38j . Fitting errors are estimated using 
the jackknife method. 

To increase statistics, we do multiple measurements on 
each configuration. For each measurement, the source 
position ti is chosen randomly, with t 2 determined by 
t 2 = ti + At, where At is the wall-source separation 
listed in Table EYi In addition, we use different random 
numbers for the wall sources for each measurement. The 
number of measurements for each gauge configuration is 
listed in Table HI 

To study auto-correlations we bin adjacent lattices in 
the Markov chain and study the dependence of the nom¬ 
inal statistical error on bin size. Examples of the results 
are shown in Fig.[^ The notation for the operators used 
in this figure is explained in Refs. [3H1 ED] ■ We find that 
the auto-correlations increase as the lattice spacing de¬ 
creases. As one can see from Fig. the auto-correlation 
effect is about 100% for the MILC superfine lattice SI, 
while it is about 25% for the MILC fine lattice FI. In 
order to greatly reduce the effects of auto-correlations, 
we use bins of size 5 throughout our analysis. 


IV. CHIRAL EXTRAPOLATION 

Our analysis follows the same steps as in Refs. [3 [23]. 
The first step is the chiral extrapolation of the valence 
quark masses to their physical values. We extrapolate 
rux to the for fixed niy using a fitting form based 

on SU(2) SChPT, and then extrapolate niy to 
For SU(2) ChPT to be valid, we require that rrix 
niy. Hence, from the 10 valence quark masses listed in 
Eq. Q, we take the lightest four for nix (e.g. nix = 
{0.003, 0.006, 0.009, 0.012} on the fine ensemble) and 
heaviest three for niy (e.g. niy = {0.024, 0.027, 0.03} on 
the fine ensemble). In this way we satisfy lUx < niy 12. 

We begin by considering the extrapolation in nix, 
which we call the “X fit”. The actual extrapolation is 
done in Ap = p, which is the squared mass of 
the XX valence pion with taste ^5 (i.e. the Goldstone 
pion). Eor the physical value of this quantity we take 
Xp = = (0.158 GeV)2. At next-to- 

leading order (NLO) in SU(2) SGhPT, the light valence 
quark mass dependence of the B-parameters has been 


worked out in Ref. [55], and is 

R,(NLO) = ciFo(f) + C2A, (9) 

Xp 

where X = with A^ = 1 GeV, the Cj are coefficients 
to be determined, and 


Fo{i) = 1 ± 


327r2/2 


l{Xi) + {Li - Xi)i{Xi) 


( 10 ) 


is the chiral logarithm. Here Xp (Lp) is the squared 
mass of the taste B, flavor non-singlet, pion composed 
of two light valence (sea) quarks: Xp = ni^,^ p {Lp = 
mfi g). The functions i{X) and ((X) are chiral loga¬ 
rithms defined, for example, in Ref. [55] • In Eq. (101, the 
plus sign applies for i = K, 2, 3, and the minus sign for 
i = 4,5. 


The NLO fitting function is not accurate enough to 
describe the precise and highly correlated data. Hence, 
as in all our recent analyses [5113 Eg, we add higher order 
terms to the fitting function: 


Hj (NNNLO) = ciEo(j) -f C2A -h -h CiX"^ In^(X) 

-|-CsA^ ln(A)-|-cgA^. (H) 


The three terms X^,X^\n'^{X) and A2ln(A) are the 
generic NNLO terms in continuum chiral perturbation 
theory. We also add a single analytic NNNLO term pro¬ 
portional to X^. We use a similar fitting function for 
the X fits of gold-plated combinations, except that, by 
construction, there are no NLO chiral logarithms: 

G*(NNNLO) = Cl + C 2 X + csX^ + ciX^ ln^(A:) 

+ C5X‘^ln{X)+ceX\ ( 12 ) 


We have found that adding yet higher order terms in the 
chiral expansion does not improve the fits to either the 
Bi or Gp 

Since we have only four data points for the X fit, we 
use the Bayesian method |39j . and place constraints on 
the three higher-order fitting parameters C 4 _ 6 . Our 
prior information is that these coefficients are of 0{1). 
We thus first impose the constraints 04-5 = 0 ± 1. If the 
resulting fits have y^/d.o.f. < 1, then we accept them. If 
not, we try the less restrictive constraints C 4_6 = 0 ± 2. 
Again, we accept hts with y^/d.o.f. < 1, but otherwise 
fit again using C 4_6 = 0 ± 4. In all cases this leads to fits 
having y^/d.o.f. < 1. In this discussion, the x2 that is 
minimized is the augmented version: 


= + Xnr 


2 _ (G ~ 

Xprior ~ 2 -^ ,j2 

i=4 * 


(13) 

(14) 
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(a) FI (b) SI (c) U1 


FIG. 1. G23 (evaluated at renormalization scale — 1/a as a function of T = t — ti. Green diamonds show on the FI ensemble 
with {amx,amy) = (0.003,0.03). Blue pentagons are results from the SI ensemble with {amx,amy) = (0.0018,0.018). Brown 
squares are results on the U1 ensemble with {amx,amy) = (0.0014,0.014). The fit ranges (and resulting central values and 
error bands) are shown by the horizontal lines. 



FIG. 2. G 45 as a function of T = t — t\ at p = 1/a. The 

convention for symbols is as in Fig. 


where we set = 0 and at = 1, 2,4. These fits are done 
using the full correlation matrix, and have acceptable 
values of x^- 

Having determined the parameters ci_6, we extrapo¬ 
late the results to the physical point vrix = , and si¬ 

multaneously remove (by hand) the lattice artifacts that 
lead to taste symmetry breaking in pion masses. Specif¬ 
ically, within the chiral logarithm To(*) ’"^6 set Xb and 
Li to their physical values, as explained in Ref. [7]. In 
this way we are using knowledge from SChPT to remove 
a significant source of discretization errors. Note that 
this correction applies to the Bi but not to the G^-, since 
the gold-plated combinations have no chiral logarithms 


at NLO. 

Examples of the X fits are shown in Figs, ^[a) ^_ 

and f ’a) for Bk, G23 and G45, respectively. In all these 
fits it was sufficient to use the narrowest range of the 
Bayesian priors {(Ji = 1) in order to obtain good fits. 
We note that the statistical errors appear larger in the 
results for G 23 because of the finer vertical scale. The 
figures emphasize the fact that the extrapolation in Xp 
is relatively short. Thus the dependence on SChPT is 
relatively mild, except for the taste-breaking correction 
that we make to Bk- 

In order to estimate the systematic uncertainty in the 
X fits we consider two variations in the fitting scheme. 
The first error is obtained from the changes in the Bi and 
Gj when the prior widths a a are doubled. The second 
is obtained by repeating the fits keeping only one NNLO 
term. 


Rx(NNLO) =ciEo(R:)-bC2X-bC3X2, (15) 

Gi(NNLO) =ci -bc 2 X-bc 3 X 2 , (16) 

and using the eigenmode shift (ES) method introduced 
in Ref. [38] . The ES method tunes the fitting function in 
the direction of the eigenvectors of the covariance matrix 
corresponding to the small eigenvalues, with small shift¬ 
ing parameters rj that are constrained by the Bayesian 
prior condition: ?; = 0 ± cr^. We set cr,, from the size of 
the neglected highest order term in the fitting function, 

cr^ = 0.006 « X2(ln(X))^ (17) 

where X « 0.02. 

The total systematic error from the X fits is then ob¬ 
tained by adding these two error estimates in quadrature. 
The resulting errors are discussed in Sec. jVTj 

We next extrapolate rUy to using the three heav¬ 

iest values of the valence quark masses. This we denote 
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b 

1.0 
0.8 
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2 4 6 8 10 

bin size 

(a) [P X P][P X P]^ 

FIG. 3. Statistical errors for bare three-point functions as a function of bin size. The operators are (a) COpi*, and (b) 0^2*i 
respectively, at T = 20 (FI) and T = 30 (SI). (Red) circles are the results on FI ensemble, with {am^,amy) = (0.003,0.03); 
and (blue) crosses are the results on SI ensemble, with {arux, aruy) = (0.0018,0.018). 


FI 

SI 



S 


2 4 6 8 10 

bin size 

(b) [PxP][PxP]„ 





Yp (GeV^) 


(b) Y-fit for Bk 


FIG. 4. (a) X fits and (b) Y fits for Bk evaluated at /r = 1/a on the FI, SI and U1 ensembles. The valence strange-quark 

masses are amy = 0.03, 0.018 and 0.014, respectively. Lattice results are shown with circles (green, blue and brown for FI, SI, 
and Ul, respectively) and are ordered vertically as shown in the legend. Extrapolated results are shown with [green] triangles 
(FI), [blue] diamonds (SI) and [brown] pentagons (Ul). For the X fit, the extrapolated results lie below the curves because of 
the removal of taste-breaking effects, as described in the text. 


the “Y fit”. We expect the Bi and Gj to be smooth, an¬ 
alytic functions of Yp, since the strange quark is far from 
the chiral limit. Empirically, linear fitting works very 
well, as illustrated in Figs. m |l^[b)| and |(]|[b)| To avoid 
the problem of small eigenvalues, we use uncorrelated fit¬ 
ting for the Y fits. In all cases, fits are stable and the fit 
parameters are consistent across all lattices with a given 
nominal lattice spacing (within the statistical uncertain¬ 
ties). To estimate the systematic error in the results of 
the Y fits, we repeat the fits using a quadratic function 
of Yp. The changes in the final results for Bi and Gj are 
then taken as the systematic error. 


V. CONTINUUM-CHIRAL EXTRAPOLATION 


The outputs of the extrapolations in valence masses 
are values for the ^-parameters and gold-plated combi¬ 
nations on each ensemble, for continuum operators eval¬ 
uated at the renormalization scale p = 1/a. In order to 
compare these results and extrapolate them to the con¬ 
tinuum limit, and to physical sea-quark masses, we must 
use renormalization group (RG) evolution to evolve to a 
common scale. The standard choices for this scale in the 
literature are p = 2 GeV and p = 3 GeV, and we present 
results for both. Since we use one-loop matching, to do 














































FIG. 5 . (a) X fits and |(b)| Y fits for G23. Notation as in 
taste-breaking correction. 



Yp (GeV^) 


(b) Y-fit for G 23 

Fig. lil except that for the gold-plated combinations there is no 



(a) X-fit for G 45 


FIG. 6. 


(b) Y-fit for G 45 

|(a)|X fits and|(b)| Y fits for G45. Notation as in Fig. 


the running consistently we need the continuum two-loop 
anomalous dimension matrix. This has been calculated 
in Ref. [dl] for a particular choice of evanescent operators. 
Because of this, it is essential that our lattice-continuum 
matching uses the same set of evanescent operators, as is 
indeed the case in Ref. [30] • Some technical issues arise 
in the RG running; these are described in Ref. m along 
with our resolutions. 


We present our results for Bk and the gold-plated 
combinations Gi at the two renormalization scales in Ta¬ 
bles |V] and VI Statistical errors range from the percent 
level to an order of magnitude smaller. We have also ob¬ 
tained results for the Bj {j = 2 — 5) but do not show 
these as they are not used in our final analysis. 

The final step of our analysis is to do a simultane¬ 
ous extrapolation to the physical values of the sea-quark 
masses and to the continuum limit. We call this proce¬ 


dure “the continuum-chiral extrapolation”, although this 
name is slightly misleading as the valence chiral extrap¬ 
olation has already been done. As substitutes for sea 
quark masses, we use Lp and Sp, which are, respec¬ 
tively, the squared masses of taste -^5 (Goldstone) pions 
composed of two light sea quarks {ll) and two strange 
sea quarks (ss). They are extrapolated to their physical 
values, which we take to be = (0.1349766 GeV)^ for 
Lp and = (0.6858GeV)^ for Sp [37] . 


We expect the dependence of the Bi and Gj on Lp, 
Sp and to be analytic, with terms organized accord¬ 
ing to standard SChPT power counting. At NLO, the 
only term in SChPT that could violate this expectation 
is the chiral logarithm. This is absent for the Gj. For 
the Bi, as shown by Eq. (10), the only logarithms that 
appear have the schematic dependence {Lp + o?) log 
and Xp logXp. Since Xp is set by hand to its physical 
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TABLE V. Bk and gold-plated combinations for — 2 GeV on each lattice listed in Table|I] The superscripts indicate whether 
broadened Bayesian priors have been used in the X-fits: f implying C 4_6 = 0±2, while J implying C 4-6 = 0±4. Results without 
superscripts are obtained with C4_6 = 0 ± 1. 


ID 

Bk 

G 21 

G 23 

G 24 

G 45 

Cl 

0.5484(55) 

0.995(ll)f 

1.4140(10) 

0.6205(19) 

1.1836(7) 

C2 

0.5528(56) 

0.993(11)''^ 

1.4119(11) 

0.6232(22) 

1.1824(8) 

C3 

0.5673(52) 

0.975(10)* 

1.4098(9) 

0.6256(20) 

1.1819(7) 

C3-2 

0.5715(51) 

0.974(9)* 

1.4105(9) 

0.6229(16) 

1.1823(6) 

C4 

0.5641(54) 

0.987(11)* 

1.4113(9) 

0.6291(19) 

1.1822(6) 

C5 

0.5677(46) 

0.976(8)* 

1.4082(8) 

0.6264(15) 

1.1834(5) 

FI 

0.5294(43) 

1.0451(69) 

1.4000(10) 

0.6092(19) 

1.2003(9) 

F2 

0.5451(35) 

1.0281(59) 

1.3985(6) 

0.6088(11) 

1.1993(5) 

F3 

0.5226(49) 

1.053(10)* 

1.4015(11) 

0.6119(21) 

1.1991(10) 

F4 

0.5255(30) 

1.0366(64)* 

1.4033(7) 

0.6050(12) 

1.2008(6) 

F5 

0.5388(43) 

1.0322(84)* 

1.3995(9) 

0.6101(17) 

1.1997(8) 

F6 

0.5472(59)^ 

1.014(11)* 

1.3991(13) 

0.6123(23) 

1.2018(9) 

F7 

0.5392(35) 

1.0394(59) 

1.3953(7) 

0.6130(12) 

1.1992(6) 

F9 

0.5501(16) 

1.0258(30) 

1.3976(3) 

0.6093(6) 

1.1991(3) 

SI 

0.5359(38) 

1.0531(55) 

1.4140(10) 

0.5858(17) 

1.2288(8) 

S2 

0.5361(36) 

1.0423(57) 

1.4116(9) 

0.5833(13) 

1.2278(8) 

S3 

0.5261(41)t 

1.0625(79)* 

1.4184(13) 

0.5842(20) 

1.2278(10) 

S4 

0.5204(33) 

1.0621(60) 

1.4124(10) 

0.5820(18) 

1.2277(8) 

S5 

0.5384(36) 

1.0446(55) 

1.4110(8) 

0.5835(12) 

1.2287(8) 

U1 

0.5325(70) 

1.056(11) 

1.4302(28) 

0.5718(39) 

1.2539(19) 


TABLE VI. Bk and gold-plated combinations for fi — 3 GeV on each lattice listed in Table [I] The convention for f and J is the 
same as Table lYl 


ID 

Bk 

G 21 

G 23 

G 24 

G 45 

Cl 

0.5298(53) 

0.951(10)* 

1.3942(8) 

0.5713(18) 

1.1468(6) 

C2 

0.5341(54) 

0.950(10)* 

1.3926(8) 

0.5738(20) 

1.1459(6) 

C3 

0.5481(50) 

0.9323(97)* 

1.3911(7) 

0.5760(19) 

1.1455(5) 

C3-2 

0.5521(49) 

0.9308(89)* 

1.3915(7) 

0.5735(14) 

1.1458(5) 

C4 

0.5449(52) 

0.944(10)* 

1.3921(7) 

0.5792(18) 

1.1457(5) 

C5 

0.5484(45) 

0.9327(80)* 

1.3898(6) 

0.5767(14) 

1.1467(4) 

FI 

0.5115(42) 

0.9991(66) 

1.3829(7) 

0.5610(18) 

1.1594(7) 

F2 

0.5266(34) 

0.9828(56) 

1.3817(5) 

0.5606(10) 

1.1586(4) 

F3 

0.5049(47) 

1.0069(96)* 

1.3840(9) 

0.5634(19) 

1.1584(8) 

F4 

0.5077(29) 

0.9909(61)* 

1.3853(5) 

0.5571(11) 

1.1597(4) 

F5 

0.5205(42) 

0.9867(80)* 

1.3825(6) 

0.5618(16) 

1.1589(6) 

F6 

0.5287(57)* 

0.969(10)* 

1.3821(10) 

0.5639(21) 

1.1604(7) 

F7 

0.5210(34) 

0.9935(56) 

1.3793(5) 

0.5645(11) 

1.1585(5) 

F9 

0.5314(16) 

0.9806(29) 

1.3811(3) 

0.5611(5) 

1.1585(2) 

SI 

0.5178(37) 

1.0068(53) 

1.3927(8) 

0.5394(16) 

1.1806(6) 

S2 

0.5179(34) 

0.9965(55) 

1.3909(7) 

0.5372(12) 

1.1798(6) 

S3 

0.5083(39)* 

1.0158(76)* 

1.3960(10) 

0.5380(19) 

1.1798(8) 

S4 

0.5028(31) 

1.0153(58) 

1.3916(8) 

0.5359(16) 

1.1797(6) 

S5 

0.5202(34) 

0.9986(53) 

1.3905(6) 

0.5373(12) 

1.1805(6) 

U1 

0.5145(67) 

1.009(10) 

1.4044(21) 

0.5266(36) 

1.1991(14) 


value, the dependence it contains is removed. The 
remaining dependence on Lp and is analytic, and in 
fact also is removed by hand when we set Lp to its phys¬ 
ical value and to zero. Chiral logarithms of higher 
order can lead to non-analyticities, or large derivatives, 
but these are numerically suppressed. Thus, to good ap¬ 
proximation, we expect all the quantities we calculate to 


be described by 

r c _ j\^2 

=di+d2^+d3^-r/^^+d4(aAQ)2. (18) 

Here Ag = 0.3 GeV and A^ = 1 GeV, and we have 
chosen to expand the term about the physical ss mass. 
When we fit our results to this form, we impose 
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Bayesian constraints on the linear terms to enforce the 
expected power counting: d 2 -i = 0 ± 2. We have also 
tried fits with broader contraints, d 2 -i = 0 ± 4, but find 
that these do not significantly change the or the re¬ 
sulting fit parameters. We find, as was the case in our 
earlier work [H El |23] that we cannot obtain a good de¬ 
scription if we include the coarse lattices. Thus we fit all 
the fine, superfine and ultrafine lattice data to Eq. (18 1 . 
We call this the Fi fit, since it is a small variation from 
the fitting function Fg in our previous work [9] (differing 
only in the offset in the d^ term). Since the number of 
configurations differ on each ensemble, errors on the fit 
parameters are obtained using a variant of the bootstrap 
method. Note that for this fit there are no correlations 
between th e dif ferent ensembles. 

In Table VII we show the results of the Fi fits to Bk 
and the Gi (renormalized at /x = 2 GeV). Plots of the fits 
are shown in Figs.[i|^ ^[a) l(l[a)[ and l] Ka)| The 

fits are qualitatively similar and of comparable quality if 
the operators are renormalized at /x = 3 GeV. To inter¬ 
pret these plots the following must be kept in mind. For 
each nominal value of a (e.g. for the fine lattices) there 
is a variation in the actual values of a and in the val¬ 
ues of Sp. This is most significant for the ensembles F 6 
and F7, which have a substantially lower strange quark 
mass than the other fine ensembles. These variations are 
accounted for in the fit (with F 6 and F7 providing a sig¬ 
nificant lever arm to determine da), but do not show up 
in these two-dimensional plots. Indeed, the points from 
F 6 and F7 are not included in the plots, while the fit 
lines for the fine and superfine ensembles are shown with 
a and Sp set to their average values (excluding ensem¬ 
bles F 6 and F7 for the fine lattices). Thus, even if the fit 
were perfect, the fit lines would not pass through any of 
the points, except for the ultrafine case. Because of this, 
the fits appear slightly worse than they actually are; the 
real indicator of goodness of each fit is the quoted value 
of x^/d.o.f.. 


The fits indicate that the dependence on the strange 
sea-quark mass is very weak for all five quantities, with 
I da I ^ 1. For the gold-plated combinations, the depen¬ 
dence on the light sea-quark mass is also weak, much 
weaker within our range of parameters than the depen¬ 
dence on a. Only for Bk does the variation with Lp 
have a similar magnitude to that with a. The values 
of the coefficient, d 4 , indicate that scale describing 
aS effects ranges from ^ 0.3 GeV (|d 4 | ~ 1) up to 
~ 0.55 GeV (|d 4 | ^ 3.5). These scales are not unusual for 
discretization errors with improved staggered fermions. 
The Y^/d.o.f. of these fits is reasonable for Bk, G 21 , and 
G24n Hence, we choose the results from the Fi fits for 


^ Here we consider a value up to ~ 1.5 to be reasonable, due to 
residual correlations between configurations. We work with a bin 
size of 5, and Fig.[^shows that this can lead to an underestimate 
of the error by ~ 25% on some configurations. Consequently the 
will be overestimated by ~ 1.25^. 


our central values for these quantities. However, we can¬ 
not use the Fi results for G23 and G45, since the fit qual¬ 
ity is too poor. This is primarily due to the difficulty 
that the fits have in reproducing the dependence on a. 

To obtain reasonable fits for G23 and G45, we add 
higher order terms to the fitting function, denoting the 
new form F 4 : 


I J / A 


iV = Fi-I-fi5(aAQ) —-|-d6(aAQ) [ 


A2 

X 


-I- d7{ah.QYas -f ds,al -I- dgiaAq^ 


(19) 


where as = °(l/a). In other words, we add a subset 

of the analytic terms quadratic in Lp, Sp and a^, as 
well as two terms that include logarithmic dependence 
on a. The dy term would be the dominant source of 
a dependence were the action and operators tree-level 
G(a^) improved. In fact, our valence fermion action and 
operators are not tree-level improved, so we must include 
the pure aS d^ term as well. Nevertheless, we expect 
the tree-level contributions proportional to alone to 
be small, due to the use of HYP-smeared gauge fields. 
The ds term arises because our lattice operators are only 
matched to the continuum operators at one-loop order, 
leaving a two-loop residual discrepancy. In the F 4 fits 
we constrain d 2_9 using the Bayesian method, choosing 
the prior conditions d^-Q = 0 ± 2. Again we find that 
broadening the priors does not significantly improve the 
fits. 

The results for the F 4 fits are shown (for /x = 2 GeV) in 
Table VHI and Figs, "t 'b) f [t ^ and [Tifb)~ 


With the additional terms, we obtain reasonable values 
of x^/d.o.f for G23 and G45, and we take the resulting ex¬ 
trapolated values as our final results for these two quan¬ 
tities. For the other quantities, the fit quality is only 
slightly improved. 

As is apparent, particularly from Figs. 


If 


b) 1C 'b) and 


b) the change from Fi to F 4 fits has a very significant 


impact on the continuum extrapolation. This is primarily 
due to the dgol term, which has a rapid dependence on 
a as a —>■ 0. We note that the coefficients of this term 
in the fits to G23 and G45 are relatively large [although 
still of G(l)], and this is what leads to the large change 
in the extrapolated value between the fits. We do not 
find the F 4 fits to provide a convincing description of the 
a dependence, particularly as they depend very strongly 
on the result from the single ultrafine lattice. However, 
we think that the conservative choice is to use the better 
fit for the central value, and then to take the difference 
between the two fits as an estimate of the systematic 
error in the continuum-chiral extrapolation. The final 
results from the two fits, and the resulting estimate of 
the systematic error, are collected in Table m For G23, 
G24 and G45 this source of error dominates all others, as 
discussed in the following section. 

We have also used fit functions with additional higher- 
order terms. These do lead to mild reductions in the 
values of y^/d.o.f, but do not lead to significant changes 
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TABLE VII. Results of Fi fits to Bk and the gold-plated combinations. The renormalization scale is — 2 GeV. 



Bk 

G 21 

G 23 

G 24 

G 45 

di 

0.5390(37) 

1.0568(62) 

1.4248(10) 

0.5590(15) 

1.2567(8) 

d2 

-0.127(14) 

0.095(27) 

0.0275^3) 

-0.0097(56) 

0.0041(26) 

ds 

0.006(15) 

-0.014(25) 

0.0145^0) 

-0.0207(53) 

0.0026(25) 

d4 

0.78(25) 

-1.92^1) 

-1.799(64) 

3.21(10) 

-3.529(53) 

Bk or Gi 

0.5366(36) 

1.0585(59) 

1.4253(10) 

0.5589(15) 

1.2568(8) 

xVdof 

1.53 

1.30 

2.01 

1.08 

4.07 




(a) Fi: x^/d.o.f = 1.53 (b) F 4 : x^/d.o.f = 1.52 

FIG. 7. Gontinuum-chiral extrapolation for Bk renormalized at /r = 2 GeV. Results from the fine, superfine and ultrafine 
lattices are shown with (green) triangles, (blue) diamonds and the (brown) pentagon, respectively, (a) Fi fit; (b) F 4 fit. The 
(red) circle gives the extrapolated result. Due to the variations in values of Sp and a, the curves should not pass precisely 
through all the points. For more discussion, see text. 


in the central values compared to the F 4 fits. Thus they 
do not significantly change our estimates of systematic 
errors. For the sake of brevity, we do not display the 
results of these more elaborate fits. 

We close this section by returning to the option of di¬ 
rectly fitting the BSM i?-parameters rather than using 
the gold-plated combinations. In all cases we find that 
direct continuum-chiral fits have values of y^/d.o.f. in 
the range 3 — 5, both for Fi and F 4 (and more elabo¬ 
rate) fits. We do not fully understand this failure of the 
continuum-chiral fits, but suspect that it is related to er¬ 
rors in valence chiral extrapolation (X fits). The X fits 
are better controlled with the gold-plated combinations. 


VI. FINAL RESULTS AND ERROR BUDGET 


XI while the final error budget is given in Table XII 


As can be seen from Table |X^ the statistical errors in 
Bk and the Gi are small, ranging from ^ 0.25% to ^ 1%. 
The largest are those in G23 and G45, resulting from the 
use of the F 4 fits for the continuum-chiral extrapolation. 
We propagate the statistical errors into the using 
the bootstrap method. The larger errors in G23 and G45 
then lead to B^ and B^ having the largest statistical 
errors of the BSM i3-parameters. In all cases, however, 
the statistical errors are much smaller than those from 
systematic effects. 

We now run through the systematic errors in the or¬ 
der fisted in Table ixm The dominant error is that due 
to the combined effect of using one-loop matching and 
the continuum-chiral extrapolation. We combine these 
because the F 4 fit includes the error that results from 
perturbative truncation, and indeed this is the domi¬ 
nant contribution to the systematic error estimate, as 


In this section we discuss all sources of error, and give 
our final results for the BSM i3-parameters with their 
error budget. Because we obtain using Eq. § , we 
estimate the errors in and the Gi hrst, and then prop¬ 
agate the errors to Our hnal results for the two 

standard renormalization scales are given in Tablesp^and 


^ The result quoted here for Bx(2Ge'V) is obtained by a very 
slightly different analysis method than that we used previously 
in Ref. [^. Thus the results differ slightly, although they agree 
within the (small) statistical errors, and have almost exactly the 
same total error. 


























G24 ( 2 GeV) G23 ( 2 GeV) G21 ( 2 GeV) 




(a) Fi: x^/d.o.f = 1.30 (b) F 4 : xVd-o.f = 1.23 

FIG. 8. Continuum-chiral extrapolation results for G 21 at /i = 2 GeV. The notation is as in Fig. 




(a) Fi: x^/d.o.f = 2.01 (b) Kj: x^/d.o.f = 1.33 


FIG. 9. Continuum-chiral extrapolation results for G 23 at /r = 2 GeV. The notation is as in Fig. 
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FIG. 10. Continuum-chiral extrapolation results for G 24 at /r = 2 GeV. The notation is as in Fig. 
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CM 
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1.18 




Fine '—A- 
Superfine ■—$- 
Ultrafine ■—6- 


—«- 


0.1 


0.2 


0.1 


0.2 


(a) Fi: x^/d.o.f = 4.07 


(b) F 4 : x^/d.o.f = 1.39 


FIG. 11. Continuum-chiral extrapolation results for G 45 at /r = 2GeV. The notation is as in Fig. 


TABLE VIII. Fit results for BSM B-parameters and the gold-plated combinations obtained using E 4 -fit at p = 2 GeV. 



Bk 

G 21 

G 23 

G 24 

G 45 

di 

0.5308(99) 

1.080(12) 

1.488(14) 

0.523(12) 

1.378(14) 

d2 

-0.124(18) 

0.104(28) 

0.045(14) 

0.001(16) 

-0.013(12) 

ds 

0.005(15) 

-0.011(25) 

0.019(6) 

-0.021(6) 

0.012(7) 

di 

0.24(40) 

-0.33(27) 

2.22(75) 

0.84(63) 

3.70(80) 


-0.18(77) 

-0.66(48) 

-1.10(87) 

-0.7(10) 

1.16(78) 

dg 

0.10(10) 

-0.081(63) 

-0.29(34) 

0.03(22) 

-0.64(42) 

dr 

0.09(21) 

-0.10(14) 

1.20(40) 

0.33(33) 

2.06(42) 

dg 

0.22(21) 

-0.65(20) 

-1.80(37) 

0.98(31) 

-3.34(39) 

dg 

0.008(31) 

-0.008(21) 

0.183(60) 

0.036(50) 

0.322(64) 

Bk or Gi 

0.5285(98) 

1.082(12) 

1.489(13) 

0.523(12) 

1.378(14) 

xVdof 

1.52 

1.23 

1.33 

0.91 

1.39 
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TABLE IX. Results for Bk and Gi (renormalized at ^ = 
2 GeV) from continuum-chiral extrapolation using the Fi and 
Fi fits. Our choices for the final central values are (red 
and) underlined. A is the fractional systematic error in the 
continuum-chiral extrapolation, and is obtained from the dif¬ 
ference between the two fits. 



Fi 

Fi 

A(%) 

Bk 

0.5366(36) 

0.5285(98) 

1.52 

G21 

1.0585(59) 

1.082(12) 

2.18 

G23 

1.4253(10) 

1.489(13) 

4.26 

G24 

0.5589(15) 

0.523(12) 

6.36 

G45 

1.2568(8) 

1.378(14) 

8.79 


TABLE X. Final results for Bk and the BSM B-parameters 
at renormalization scales /r = 2 GeV and fj, = 3 GeV. The 
first error is statistical, the second systematic. 



H = 2 GeV 

/i = 3 GeV 

Bk 

0.537(4)(26) 

0.519(4)(26) 

B§ 

0.568(1)(25) 

0.525(1)(23) 

B§ 

0.382(4)(17) 

0.360(4)(16) 

Bf 

0.984(3)(64) 

0.981(3)(62) 

B? 

0.714(7)(71) 

0.751(7)(68) 


discussed above. However, one can also estimate the 
truncation error directly, following Ref. [5], by the size 
of a typical two-loop contribution: 


ABi ^ BiX al. ( 20 ) 

Here we use ag in the MS scheme evaluated at a scale 
l/oinin, where amin is our smallest lattice spacing. This 
leads to a 4.4% relative error. To avoid double-counting, 
we take the larger of (a) the direct estimate of two-loop 
effects (4.4%) and (b) the difference between Fi and F 4 
fits. In essence this method is using the F 4 fit to give 
an estimate of the uncertainty in the coefficient of the 
ttg term, except that we do not allow this uncertainty to 
drop below unity. 

The above description applies to quantities we calcu¬ 
late directly, namely Bk and the Gi. For the derived 
quantities B^, defined in Eq. we proceed as follows. 
We vary the fit choices (for the continuum-chiral extrap¬ 
olation) for each of the components of the B^ indepen¬ 
dently, and take the largest variation from the central 
value as the systematic error estimate. If this maximum 


TABLE XL Final results for the gold-plated combinations Gi. 
Notation as in Table m 



= 2 GeV 

fi = 3 GeV 

G 21 

1.059(6)(52) 

1.012(6)(50) 

G23 

1.489(13)(66) 

1.460(14)(65) 

G24 

0.559(1)(36) 

0.515(1)(32) 

G45 

1.378(14)(123) 

1.307(14)(107) 


value lies below 4.4%, we replace the estimate with the 
direct two-loop estimate of a 4.4% error. 

We next consider the error due to the finite volume 
(FV) of the lattice. We estimate this from the difference 
between results on the C3 and C3-2 ensembles, which 
differ only in their spatial volumes. This is not entirely 
satisfactory, since we do not use coarse lattices in our fi¬ 
nal continuum-chiral extrapolation. However, we stress 
that the dominant FV error, as estimated by SChPT, 
comes from valence pions propagating to adjacent peri¬ 
odic volumes. This is because the arguments of the chiral 
logarithms of Eq. (101 are the squared masses of valence 
pions, Xb. Since on each ensemble we are extrapolating 
to the physical valence quark masses, the dominant FV 
effects are present, even though on ensembles C3 and C3- 
2 we are far from the physical values of Lp and a. In our 
calculation of Bk, we have used the comparison of doing 
the X-fits with finite- and infinite-volume SChPT forms 
as an alternative estimate of the FV error |40]. However, 
this method is not useful for the gold-plated combina¬ 
tions, since they do not contain NLO chiral logarithms. 

Our method of estimating systematic errors arising 
from the X fits has been described in Sec. IV We repeat 
the entire analysis using different priors for the X-fits, 
and using the ES method. Each leads to a change in the 
final values of the quantities of interest. We combine the 
fractional shifts in quadrature to obtain our total system¬ 
atic error. 

Our method of estimating the systematic error aris¬ 
ing from Y fits, as noted above, is to repeat the entire 
analysis (including the continuum-chiral extrapolation) 
using quadratic, as apposed to linear, functional forms. 
This differs slightly from the estimate we used in Ref. [9] , 
where we used the shift in the quantities on a specific 
MILC ensemble. The Y fit errors turn out to be of com¬ 
parable size to those from X fits, ranging up to 2%. 

The remaining two systematic errors are very small, 
and have essentially no impact on the total error. We 
include them for completeness. The first concerns the 
value of the pion decay constant / that we use in the 
chiral logarithms of Eq. (10). At NLO we could equally 
well use the physical value fir = 130.41 MeV m or the 
value in the chiral limit, /,r ~ 124.2 MeV [32]. In prac¬ 
tice we use / = 132 MeV (close to the physical value) 
for our central value, and repeat the entire analysis us¬ 
ing / = 124.2 MeV (the chiral-limit value) to estimate 
the systematic error. In fact, only Bk is sensitive to 
this choice, since the gold-plated combinations contain 
no NLO chiral logarithms. Thus the 0.1% error that re¬ 
sults in Bk propagates unchanged into all of the BSM 
R-parameters. 

Finally, the parameter we use to set the scale, ri, has 
an error which propagates into all the final results. To es¬ 
timate this, we repeat the entire analysis with the central 
value for ri replaced by ri ±(7^, and quote the maximum 
difference in each quantity as the systematic error. The 
resulting errors are very small (~ 0.1 — 0.35%), reflecting 
the fact that the R-parameters are dimensionless. 
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TABLE XII. Error budget for the Bi and Gj evaluated at renormalization scale /r = 2 GeV. All entries are in percent. 


cause 

Bk 

G 21 

G 23 

G 24 

G 45 


B^ 

to 

Bf 

dG,SUSY 

-^3 

method 

statistics 

0.67 

0.56 

0.87 

0.27 

1.02 

0.25 

1.00 

0.27 

0.98 

0.66 

see text 

J matching 1 

1 cont-extrap. J 

4.40 

4.40 

4.40 

6.36 

8.79 

4.40 

4.45 

6.36 

9.63 

4.40 

(Ri vs. F 4 ) or al (Ul) 

finite volume 

0.73 

0.17 

0.05 

0.43 

0.04 

0.56 

0.52 

0.99 

1.02 

0.60 

(C3) vs. (C3-2) 

X-fits 

0.05 

0.40 

0.45 

0.02 

0.96 

0.36 

0.34 

0.37 

1.23 

0.60 

change Bayes, prior & fit method 

Y-fits 

2.07 

2.11 

0.32 

0.48 

1.12 

0.00 

0.32 

0.48 

1.59 

0.22 

linear vs. quad. 

fn 

0.10 

0 

0 

0 

0 

0.10 

0.10 

0.10 

0.10 

0.10 

132MeVvs. 124.2 MeV. 

ri 

0.35 

0.28 

0.11 

0.16 

0.21 

0.07 

0.17 

0.09 

0.30 

0.01 

errors due to ri ambiguity. 

Total 

4.93 

4.91 

4.44 

6.39 

8.92 

4.45 

4.51 

6.47 

9.90 

4.49 



VII. COMPARISONS AND OUTLOOK 


In Table |XIII| and Fig. we compare our results 
for the B-parameters to those from other collaborations. 
This is done at /x = 3 GeV since results from all collabora¬ 
tions are available at this choice of renormalization scale. 
The RBC-UKQCD collaboration uses Nf = 2-1-1 light 
flavors of domain wall quarks, and NPR for the match¬ 
ing between lattice and continuum theories. In 2012, 
RBC-UKQCD used the RI-MOM scheme for this match¬ 
ing [21], while the preliminary 2015 results are obtained 
using several RI-SMOM schemes, in the spirit of Ref. m- 
Both schemes are connected to the MS scheme using one- 
loop perturbation theory. The ETM collaboration uses 
twisted-mass Wilson quarks. The original results from 
2012 were with Nf = 2 light sea quarks and a quenched 
valence strange quark m, while the 2015 results are from 
an Nf = 2 -|- 1 -l- 1 simulation including both strange and 
charmed sea quarks |3]. Both ETM calculations match 
lattice and continuum operators using NPR in the RI- 
MOM scheme. 

Both RBC-UKQCD and ETM results are quoted us¬ 
ing the so-called SUSY basis of BSM four-fermion opera¬ 
tors |12|. The only BSM R-parameter which differs from 
that in the basis of Buras et al. (Ref. [3T|) that we use 
is Rs, 


dSUSY ^ oBuras , ^ oBuras 

^3 =- 2^3 + 2 '^^ 


( 21 ) 


We use this equation to determine our result for R|^®^ 
quoted in Table |XlII| 

For completeness, we note that our 2013 results for the 
BSM R-parameters (Ref. [2^ 1 are superseded and cor¬ 
rected by our present resultslj We now have significantly 
more ensembles, allowing a better controlled continuum- 
chiral extrapolation. This addition required us to change 
from Fi hts to R 4 hts for G 23 and G 45 , which, as shown 
above, signihcantly changes the central values for these 
quantities. In addition, we found an error in our RG run¬ 
ning due to the use of an incorrect two-loop contribution 


This does not apply to Bk , for which our present result is essen¬ 
tially the same as that from Ref. m- 


to the pseudoscalar anomalous dimension [needed for the 
denominators of the BSM R-parameters—see Eq. ©]■ 
Correcting this error leads to ~ 5% reductions in all the 
BSM R-parameters. A detailed description of the effect 
of these changes is given in Ref. [24j. The overall ef¬ 
fect is that our new results for R 2 , R 4 and R 5 

are reduced by about 5%, 3%, 5% and 12%, respectively, 
compared to those in Ref. [55] . 

Table |XIII| shows that the results for Bk, R 2 and 

rSUSY 

are consistent across all calculations, with all re¬ 
sults having comparable errors. By contrast, there are 
significant differences for R 4 and R 5 , as one can see 
most clearly from Eig. |12| The preliminary results from 
RBC-UKQCD (2015) using the intermediate RI-SMOM 
schemes are consistent with our results, while those us¬ 
ing the intermediate RI-MOM scheme [RBC-UK (2012), 
ETM (2012) and ETM (2015)] differ significantly. For 
example, the ETM (2015) results for R 4 and R 5 differ 
from our results by 2.6a and 3.2cr, respectively. 

The pattern of results in the Table suggests that the ul¬ 
timate source of these differences may well be the match¬ 
ing from lattice matrix elements to those in the contin¬ 
uum MS scheme. In our calculation, this error is due 
to the truncation of matching factors at one-loop order. 
For R 4 and R 5 (the two R-parameters which differ from 
the results obtained using the RI-MOM scheme) our er¬ 
ror estimate is taken as the difference between fits using 
Fi and R 4 fit forms (see Figs. 10 and 11). While we 


consider this to be a conservative estimate, we cannot 
rule out that it is an underestimate due to unexpect¬ 
edly large terms in the matching factors. In the case 
of the calculations using the NPR method, the signifi¬ 
cant differences between results obtained using RI-MOM 
and RI-SMOM schemes indicate an underestimate of the 
associated systematic errors. This could be a problem 
specifically related to the RI-MOM scheme, where one 
must subtract unwanted contributions from pion poles, 
a source of systematic errors absent in the RI-SMOM 
schemes [35]. Or it could be due to large truncation er¬ 
rors in the relation between one or both of these schemes 
and the MS scheme. 


Clearly these issues require further investigation. One 
possibility is for all the calculations to use the same in¬ 
termediate scheme such as RI-SMOM and then to di- 
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TABLE XIII. Comparison of the BSM B-parameters at renormalization scale /r = 3 GeV obtained using different fermion 
discretizations. The RBC-UKQCD results using domain-wall fermions are RBC-UK (2012) [21] and the preliminary results 
(with incomplete error budget) of RBC-UK (2015) [25]. The ETM collaboration results using twisted-mass fermions are ETM 
(2012) [13j and ETM (2015) |3]. N.A. means “not available”. 



SWME (this work) 

RBC-UK (2012) 

RBC-UK (2015) 

ETM (2012) 

ETM (2015) 

Bk 

0.519(4)(26) 

0.53(2) 

0.53(1) 

0.51(2) 

0.51(2) 

B2 

0.525(1) (23) 

0.43(5) 

0.49(2) 

0.47(2) 

0.46(3) 

dSUSY 

-03 

0.773(6)(35) 

0.75(9) 

0.74(7) 

0.78(4) 

0.79(5) 

oBuras 

-O3 

0.360(4)(16) 

N.A. 

N.A. 

N.A. 

N.A. 

B4 

0.981(3)(62) 

0.69(7) 

0.92(2) 

0.75(3) 

0.78(5) 

B5 

0.751(7)(68) 

0.47(6) 

0.71(4) 

0.60(3) 

0.49(4) 



■—e—' 

SWME (2015) 


■—e— 

SWME (2014) 


■e- 

RBC-UK (2015, RI-SMOM) 

■—e—■ 


RBC-UK (2012, RI-MOM) 

— e— 


ETM (2015, RI-MOM) 

-e- 


ETM (2012, RI-MOM) 


0.6 0.8 1 1.2 1.4 1.6 

B4 

(a) B4(3GeV) 


■ -e-■ SWME(2015) 

■ -e-■ SWME(2014) 

-e- RBC-UK (2015. RI-SMOM) 

■ —e —■ RBC-UK (2012, RI-MOM) 

^-e-^ ETM (2015, RI-MOM) 

-e- ETM (2012, RI-MOM) 

0.4 0.6 0.8 1 1.2 1.4 

B5 

(b) B5(3GeV) 


FIG. 12. Comparison of results for B 4 and B 5 at /r = 3 GeV. The references for the points are, proceeding from top to bottom, 
this work (SWME 2015), [24] (SWME 2014), [25] (RBC-UK 2015), [H] (RBC-UK 2012), [3] (ETM 2015) and [l3] (ETM 2012). 


rectly compare results in that scheme. This reduces 
the dependence on perturbation theory as one does not 
need to to match to the MS scheme. One would still 
need to evolve between different scales in the RI-SMOM 
scheme, but this could also, ultimately, be done non- 
perturbatively [U]. To these ends we are pursuing the 
implementation of NPR using staggered fermions gSHlZI- 
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